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There is an abundance of empirical evidence in the numerical relativity literature that the form in 
which the Einstein evolution equations are written plays a significant role in the lifetime of numerical 
simulations. This paper attempts to present a consistent framework for modifying any system of 
evolution equations by adding terms that push the evolution toward the constraint hypersurface. 
The method is, in principle, applicable to any system of partial differential equations which can be 
divided into evolution equations and constraints, although it is only demonstrated here through an 
application to the Maxwell equations. 
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I. INTRODUCTION 

The data analysis needs of the LIGO and LISA exper- 
iments have driven an interest in finding numerical solu- 
tions to the Einstein equations in astrophysically relevant 
situations. Wave forms from such numerical simulations 
can, once generated, be used as templates for finding 
gravitational wave signals in the experimental noise. 

In order to find numerical solutions, however, one must 
make a 3+1 space-time split of the Einstein equations fol- 
lowing some variation of the ADM formalism pj . When 
using free evolution techniques, which is the most com- 
mon practice in numerical relativity today, a plethora of 
authors have found that the exact form of the evolution 
equations, in particular how the constraints are substi- 
tuted into the evolution equations, plays a critical role in 
how stable the evolution will be (e.g. |1 H S S H ) . Al- 
though making such modifications to the evolution equa- 
tions remains more of an art than a science, some authors 
have proposed well defined methods for making stability- 
friendly modifications. 

Historically, Detweiler Q seems to be one of the first 
authors to consider explicitly a method for adding terms 
proportional to the constraints to the evolution equations 
with the goal of improving simulation stability. He looked 
at the evolution of the constraints, and tried to add terms 
to the evolution equations of the fundamental variables 
that would drive the constraints to zero as the evolution 
progressed. This approach is very similar to the approach 
presented here. 

Where the approach differs, however, is significant. 
Whereas Detweiler needed to add terms such that a 
negative-definite operator appeared on the right hand 
side of the evolution equation for the constraints, the 
method proposed here generates a general term of the 
correct sign which is added to terms of unknown sign 
that come from the underlying formalism. Under special 
circumstances, Detweiler showed that his method does in 
fact generate a negative-definite operator by making use 
of a free parameter introduced in his correction terms. 
The method proposed here also uses a free parameter in 



arguing that the constraints can be evolved toward zero, 
with the difference that the argument here is more gen- 
eral but without a guarantee that the overall sign of the 
the constraint evolution has the correct sign. 

The A-system approach, proposed by Brodbeck, Frit- 
telli, Hiibner, and Reula, embeds the Einstein equations 
into a larger system of symmetric hyperbolic equations. 
The Einstein equations are coupled to the larger sys- 
tem by adding terms to the Einstein evolution equations 
which are zero on the constraint hypersurface, with the 
extra degrees of freedom defined such that they are for- 
ever zero if the Einstein constraints and constraints on 
the new degrees of freedom are satisfied in the initial 
data. A linear order analysis of this system suggests that 
for "small enough" deviations from the constraint hyper- 
surface, the system should asymptote to a solution of the 
Einstein equations 

Finally, Yoneda and Shinkai [9| performed detailed 
eigenvalue calculations to determine the theoretical signs 
of the constraint evolution when linearizing around a 
Minkowski background spacetime. While their program 
was intensive and thorough, their results depend on the 
validity of a perturbation expansion around a fixed space- 
time. In that limit, they evaluate the sign of the eigenval- 
ues of their expansion matrix, and use those signs to pre- 
dict the goodness of correction terms which are strictly 
linear in the constraints and spatial derivatives of the 
constraints. Although one of their more promising terms 
was used to good effect in a numerical experiment [Toj . 
their analysis remains strictly valid only in the pertur- 
bation regime around Minkowski space, and, at present, 
they have only examined terms linear in the constraints 
and derivatives of the constraints. 

The method here, in contrast, attempts to make the 
constraint surface an attractor for the evolution equa- 
tions without appealing to perturbation theory in any 
way and without extending the system of equations. The 
results will show it successful in reaching this goal for the 
simple case of the Maxwell equations, but cannot confirm 
with certainty that the results will generalize to non- 
linear general relativity. The primary problem is that 
the correction terms proposed here have the form 
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dt 



constraint 2 = (orig rhs) — (correction) 



with a correction term of definite sign subtracted from a 
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term of indefinite sign coming from the underlying the- 
ory. The success of the method in the case of the linear 
Maxwell equations makes a similar result in the linearized 
Einstein equations likely but does not guarantee success 
for full, non-linear general relativity. This will be dis- 
cussed more completely in Section llVI 

In this paper, following earlier work by Knapp, Walker, 
and Baumgarte (KWB) |ll|, I will apply the method 
to Maxwell's equations as a proof-of-concept test. The 
work by KWB showed that some qualitative features of 
two popular ADM formulations of general relativity are 
manifest in the simpler framework of electromagnetism, 
making this a natural place to demonstrate the feasibility 
of the new methods presented here. 
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to the equations of motion, then 
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F(x,v) - K 



(dE 2 



V dx 



dE 2 
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(7) 



shows how the constraint evolves under the new equa- 
tions of motion. For large K, the term in brackets 
will dominate the evolution of the constraint, constantly 
pushing it toward its minimum value of zero. 

Numerical results, integrating forward in time with the 
Euler method, are in agreement with this analytic result. 



II. THE THEORY 



B. PDEs 



I will introduce the idea here by first considering the 
ordinary differential equation example of the simple har- 
monic oscillator, and then considering the more interest- 
ing case of partial differential equations. 



A. ODE: Harmonic Oscillator 

The equations of motion for a harmonic oscillator (with 
uj = 11 are 



dx 

~dt ~ V 
dv 

~r = — x 
dt 

and the constraint is the shifted energy function 



E(t) = x 2 + v 2 - E 



(1) 
(2) 

(3) 



with E = x 2 (0) + v 2 (0), which should be zero- valued 
at all times. In a numerical evolution, one would like 
to evolve the fundamental quantities x and v, and, after 
doing so, calculate the constraint E(t) from x and v. 
Since E is derived, on the other hand, it makes sense to 
look at its time evolution in terms of the time evolution 
of x and v using the chain rule. Furthermore, E(t) is not 
bounded on either side, so consider instead E 2 (t), which 
should also always be zero. Applying the chain rule gives 
the elementary result that 



dE 2 
~~dt 



F(x,v) 



dE 2 dx dE 2 dv 



dx dt dv dt 



(4) 



Now we cannot change the two partial derivative factors 
since they are fixed by the form of the energy expression, 
but, as noted before, we are free to change the equations 
of motion without changing the physics by adding multi- 
ples of the constraints. Let K be a positive constant. If 
we make the changes 



dx 
It 



v - K 



dE 2 
dx 



(5) 



Having seen the method in the simple case of ODEs, 
consider now the the more interesting case of PDEs with 
one constraint equation. Let s be the state vector for a 
system. Define 



S n (t,x) 



ds r . 



dt 



(8) 



to be the right hand side of the evolution equation in 
the unmodified formalism. A general constraint C will 
depend on s and its spatial derivatives. Furthermore, the 
constraint should be satisfied at every point in space. 

These considerations motivate looking at the inte- 
grated, squared constraint C 2 = J C 2 d N x, and, instead 
of taking partial derivatives with respect to the fields, we 
need to take variational derivatives of the integrated con- 
straint when considering the analogies of Equations JSHOJ 
so that the dependence of C on the spatial derivatives of 
the fields is treated properly. 

Following this prescription, the appropriate modifica- 
tion to the equation of motion for the state vector is 



ds 



dt 



S m {t,'X-) K mn [t,x) 



sc 2 



(9) 



for some positive-definite matrix-valued function K„ 
Under this change, 



dC 2 
~dT 



D[s] 




(10) 



gives the evolution of the constraint in the modified the- 
ory Here D[s] gives the functional form of the right hand 
side of the constraint's evolution equation in the unmod- 
ified theory. For the cases considered here, I chose the 
K mn diagonal and constant. 

For systems with M constraint equations, the method 
is easily modified by taking the grand constraint func- 
tional to be 



C 2 - 



w IJ {t,x)C I Cjd 



(11) 
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with any positive definite matrix wij. Like K mni the ma- 
trix Wu, can, in principle, be a function of both space 
and time, and can have off-diagonal entries. In practice, 
however, I have only used diagonal and constant matri- 
ces. Even in the diagonal case, the matrix is necessary 
because there is no a priori reason to believe that the 
constraints have the same dimensions, and it also allows 
the different constraints to be treated with different rel- 
ative importance. In addition, since there is no natural 
scale for the grand constraint, one may always set one of 
the coefficients wjj in l(TT|l to unity. 



III. APPLICATION TO MAXWELL 

This section takes Maxwell's equations as a concrete 
example of a system of partial differential equations sub- 
ject to a pointwise constraint. Following KWB [ijj, I 
consider two ways of writing Maxwell's equations for the 
vacuum in terms of the vector potential A;. The first 
system, called System I, uses the evolution equations 



d t A l = -Ei - diij) 

d t E l = -djdjAi + didjAj 



and the constraint 



C E = diEi = 0. 



(12) 
(13) 



(14) 



The second system introduces the additional field T de- 
fined by 



r = diAi 



(15) 



to eliminate mixed derivatives in (|13j) . The evolution 
equations for System II are 



BtEi = -djdjAi + diT 
d t T = -didtip 



(16) 
(17) 



and i)12fl . Both systems use a gauge consistent with 

d t tp = -d^ = -r (18) 

using the first equality for System I and the second for 
System II. 

A. System I Evolution Equations 

Having defined the systems, I would now like to calcu- 
late the terms required for applying the constraint finding 
method to System I. Here there is only one constraint, 
C = Ce = diEi, which is zero- valued. It depends only on 
the first derivatives of the electric field, therefore I need 
only to calculate 



SC 2 



■2di.CE 



(19) 



which modifies (|13fl • The new evolution equation for the 
electric field is 

d t Ei = -d k d k Ai + did k A k + 2K E d l C E (20) 

for an arbitrary positive Ke, while the other System I 
evolution equations i|12fl and (|18|l remain unchanged. 



B. System II Evolution Equations 

System II, unlike System I, has two constraints that 
should be enforced, the original constraint given by i|14|) 
plus the definition of T in (|15|) , rewritten as 



Cr 



diA - r = o 



(21) 



to make it zero-valued. This provides more freedom in 
constructing the grand constraint 



C 2 = Cf 



wC? 



(22) 



where one does not necessarily have to treat the con- 
straints on equal footing. In this case, the total con- 
straint depends additionally on T and first derivatives of 
the vector potential. In addition to l|19(l . which is still 
valid, I need 



SO 1 



= -2wdiC v 



-2wCr 



(23) 
(24) 



to enforce the definition of V. 

Applying these correction terms to the evolution equa- 
tions gives the new equations of motion 

d t Ai = -E l -d l ^ + 2wK A d l C T (25) 
dtEi = -d k d k A i + d i T + 2K E d i C E (26) 
d t T = -d k d k i> + 2wK T C r (27) 

which, combined with 118(1 . form a complete system. The 
constants Ka and Kr are arbitrary but positive. 

C. Propagation of Constraints 

In ^ij, KWB examined the evolution equation for Ce 
in both systems. They showed that for System I, the 
constraint does not evolve in time, and that in System 
II, the constraint obeys a wave equation. They did not 
have reason, however, to consider the time evolution of 
the secondary constraint Cr- I extend their results here 
by showing that the secondary constraint also satisfies 
a wave equation in the unmodified case. Furthermore, 
I demonstrate the improved behavior of the constraints 
under the modifications proposed here. 

Calculating the first time derivatives of the constraints 
is easily accomplished by taking the time derivatives of 
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(|14f> and i|21|) and replacing the time derivatives that ap- 
pear on the right hand sides of the equations by the evolu- 
tion equations l(2T)|) , lt2~5*)) , and (|2"?|) . This gives the results 



d t C E 
d t Cr 



-C E 



C r - 2K E C E ] (28) 
+ 2w [K A dfC r - K r C r ] (29) 



which can be viewed as valid for both systems if Cr is 
taken to be identically zero for System I. 

To see that KWB have a wave equation for the sec- 
ondary constraint, set all of the Ks to zero in (|28H29|) to 
eliminate the modifications, and take the time derivative 
of 123- The result 



d 2 t c r 



dfCr 



(30) 



follows immediately. 

Of greater interest here, however, is an analysis on the 
modified equations I|28l - I29l) in their first derivative form. 
The equations are linear, so they admit a Fourier analysis 
by substituting a plane wave solution e lkx into the right 
hand sides. After substituting, the resulting equations 

d t C E = k 2 C T ~2K E k 2 C E (31) 
d t C r = -C E - 2w [K A k 2 + K T ] C r (32) 

retain the terms that gave the KWB wave equations for 
the constraints, but have additional terms that look like 
they provide exponential decay. This system of equations 
is, in fact, simple enough for Maple to solve analytically 
for general values of the Ks in one dimension. The solu- 
tion, which is too long to display in detail here, consists 
of a sum of terms with the form 



exp 



-A + ± 



where 



ff{k 2 ) = {K E ± wK A )k 2 ± wK T 



(33) 



(34) 



a = ±1, and fi = /2(fc,C(0,x)) is some simple func- 
tion of k and the initial values of the constraints. Since 
fi~(k 2 ) > is manifestly positive and the radical is ei- 
ther positive or pure imaginary, the only term of this 
form that could cause anything other than exponential 

decay is —fi + \J v[(fi) 2 ~ k 2 } for parameters where 

<j[(/ 1 ~) 2 — k 2 ] > 0. Simple algebraic analysis, however, 
shows that even this term has an overall minus sign, giv- 
ing exponential decay. 

In order to make the argument more concrete, I present 
here the k = 1 solution 



C E (t,x) 
C v {t,x) 



' 3t [C E (0,x) 
~ 3t [Cr(0,x)- 



- S(x)t] 
S(x)t] 



(35) 
(36) 



for which all of the Ks are set equal to one. Here S (x) = 
C E (0, x) + Cr (0, x) is a short hand. The general solution 
can be requested from the author in the form of a Maple 
worksheet. 



Run 


Kb 


K A 


Kr 


w 


1-0 





- 


- 


- 


1-1 


i x icr 2 


- 


- 


- 


1-2 


5 x 10" 2 








II-O 











1 


II- 1 


5 x 10" 3 


5 x 10" 3 


5 x 10" 3 


1 


II-2 


i x icr 2 


1 x KT 2 


i x itr 2 


1 



TABLE I: The parameters used for the various data runs done 
on the two Maxwell systems are tabulated here. 



D. Numerical Results 

My numerical experiments on these modified systems 
of equations used an ICN integration scheme [12], and 
a Courant factor of 1/2. The spatial domain ran from 
-6 to +6, with data stored on 99 points in each coordi- 
nate direction. On the evolved fields (-Ej, Ai, ip, and 
I imposed outgoing wave boundary conditions (see 
for implementation details), and on the constraints, I im- 
posed Cj = for applicable /. All runs were performed 
on a 500 MHz Digital Personal Workstation with 1.5 GB 
of RAM. 

For each system, I ran three parameter sets, one of 
which reproduced the equations used by KWB. The full 
definitions of all of the parameter sets are found in Ta- 
ble[I] Sensible values of the parameters were easily deter- 
mined by trial and error. Choosing the values too small, 
as expected, makes little difference in the evolution, while 
choosing the values too large leads to numerical instabil- 
ities during the transient period of the evolution [T^ |. 

I followed KWB in using the analytic solution 



E4> 





8AX 2 re~ 



(37) 
(38) 



of a toroidal dipole to generate the initial data. The other 
components of the fields are zero. I chose A = A = 1, and 
the conversion from spherical to Cartesian coordinates 
was made in the code. 

The results of the System I runs are summarized in 
Figure n which shows a plot of ||Ce|| 2 vs. t. The 1-0 
(control) curve reproduces the findings of KWB that, af- 
ter an initial transient, the C E constraint does not evolve 
in time. The 1-1 and 1-2 curves, representing different val- 
ues (see Table[IJ of the parameter K E , on the other hand, 
show a modulated exponential decay. Eventually, around 
t w 200, the 1-1 case also stops decaying as rapidly, while 
the rapid exponential decay continues through the end of 
the run for 1-2. 

That the constraint in the modified case 1-1 ceases to 
evolve at some point is consistent with (|10[L which im- 
plies that the constraints will cease to evolve when the 
first term balances with the second term. From l|10fl . 
one expects that this balance will be achieved for smaller 
constraint violation when the driving parameter is larger, 
which is consistent with the results shown in Figure ^ 
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System I: ||C E || 2 vs. t 



System II: ||C E || 2 vs. t 
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FIG. 1: The h norm of the primary constraint Ce versus FIG. 2: The h norm of the primary constraint Ce versus time 
time t for three test cases. Case 1-0 has no correction terms t for three test cases. Case II-O reproduces KWB. See Tabled 



(Ke = 0). See Table [I] for the other parameter values. 



for the definitions of other parameters. 



Looking at two dimensional slices of the constraint 
data at various times, also suggests that the source of 
the modulation in the decay demonstrated by 1-1 and 
1-2 is fluctuations at the boundary, possibly caused by 
the simple boundary condition applied there on the con- 
straints. Significantly, these fluctuations are unable to 
penetrate the interior of the computational domain, un- 
like many scenarios seen in numerical relativity where 
noise from the boundary noticeably propagates inward 
from the outer boundary, eventually killing the simula- 
tion. Because I am only interested in the Maxwell equa- 
tions as a test-bed for the method, and because the mod- 
ified System I equations already perform orders of mag- 
nitude better than their unmodified counter part, I have 
not pursued this point further. 

Figure [21 shows a plot of ||Ce|| 2 vs. t for the System II 
case. Here again, the control run (11-0) reproduces the 
results of KWB, this time showing exponential decay in 
the primary constraint. Even with such an ideal result in 
the unmodified case, the modified runs II-l and II-2 show 
improvement. They represent runs with non-zero values 
(see Table HJ) of the various forcing parameters, and in 
these cases the constraint decays exponentially, but with 
a smaller characteristic time. 

Figure is a plot of ||Cr|| 2 versus time t, and shows 
similar behavior to the primary constraint in all three 
cases. The secondary constraint shows exponential decay 
in the unmodified II-O case, while showing a faster decay 
in the two modified runs, II-l and II-2. The jump in the 
graph at t = occurs because the secondary constraint 
is exactly satisfied in the initial data by construction. 



It should be noted that one likely explanation for the 
extremely favorable performance of the unmodified Sys- 
tem II is that, since the constraints satisfy a wave equa- 
tion, constraint violations propagate off of the grid. This 
is supported by the data presented by KWB, showing 
that decay rate of the constraint decreases when the outer 
boundary is moved farther out p"H. The method pre- 
sented here benefits from constraint violations propagat- 
ing off of the grid as well, but, in addition, attempts to 
damp the constraint violation throughout the grid at all 
times. 

In both systems, there is a practical limit to how large 
the forcing parameters may be chosen without creating 
an instability. I believe that its origin can be seen in 
Q33JI. which shows how the dispersion relation for the 
Fourier components of the solutions is modified by the 
new terms. This modification is effectively a modification 
of the propagation speed of those components, which, for 
a fixed Courant factor, can lead to numerical instabilities 
in the numerical scheme. When the forcing parameters 
are too large, the code crashes. The constraint blows 
up rapidly immediately before such a crash, and such 
crashes usually occur during the transient period of the 
numerical scheme. 



IV. CONCLUSIONS AND DISCUSSION 

The method proposed here does not provide a formal- 
ism for any physical system. It attempts to take any 
formalism and makes that formalism better. The analy- 



6 



System II: ||C r || 2 vs.t 




50 100 150 200 250 300 350 



FIG. 3: The h norm of the secondary constraint Cr versus 
time t for three test cases. Case II-O is the result of KWB. See 
Table^Jfor other parameter values. At t = the constraint is 
identically satisfied. 

sis does not depend on the form of the equations and the 
arguments leading to the modifications do not depend on 
the validity of any perturbative expansion. The results in 
the case of the Maxwell equations are very encouraging 
and consistent with the analytic predictions. 

Should the results derived from Maxwell's equations 



generalize to the Einstein equations, furthermore, this 
could prove extremely beneficial to the long term stability 
of simulations in numerical relativity. Fourier analysis 
of the ADM system suggests that instabilities seen in 
various 3+1 formulations are triggered by exponentially 
growing modes not seen in the initial data O n the 
other hand, at least for the case of the Maxwell equations, 
the terms added via the method proposed here lead to 
exponential damping terms in the evolution equations for 
all Fourier modes. 

It should be noted, however, that in some sense the 
Maxwell equations were an optimal case. Notice, by 
counting the number of integrations by parts required 
to calculate the correction terms to the evolution equa- 
tions, that the highest order of spatial derivatives in the 
new evolution equations will be the larger of (1) order of 
the original equations and (2) twice the order of the con- 
straints. This means that, while the order of the Maxwell 
evolution equations was not changed, the Einstein equa- 
tions will change from second order to fourth order in 
space. How this change of order effects the numeric so- 
lutions to the equations must be studied as this method 
is applied in the context of numerical relativity. 
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